function res = residual_for_agg_var(Par,vphi,mz,logagg,zH_expost,mguess)
% Description:
%   Returns the residual of labor market clearing and trade balance
%   conditions at given levels of aggregate variables, taking ex ante and 
%   ex post firm-level productivity as given
% Output:
%   res - 1 by 3 vector of residuals. 1st element is excess labor demand in
%   Home; 2nd element is excess labor demand in Foreign; 3rd element is
%   trade surplus of Home country
% Input:
%   Par - parameter structure
%   vphi - 1 by ngoods vector of parameter phi for each Home firm
%   mz - 2 by ngoods ex ante productivity matrix
%   lnagg - 1 by 3 vector of log of aggregate variables
%   zH_expost - 1 by ngoods vector of realized productivity shocks for Home firms
%   mguess - 5 by ngoods matrix of initial guesses for labor and cutoff
%       productivity. Top 4 by nfirm matrix is the guess for labor in log. 
%       Bottom row is the initial guess for cutoff productivity in log.
% Note: Aggregate variables take log form (lnagg) to accommodate searches 
% of negative numbers by fsolve

    ngoods = Par.ngoods;
    sigma = Par.sigma;
    rho = Par.rho;
    tau = Par.tau;
    kappa = Par.kappa;
    agg = exp(logagg);
    % Pre-allocate matrix to store solutions
    mlabor = zeros(4,ngoods);
    vztildeH = zeros(1,ngoods);
    
    % Specify initial guesses of variables
    mlguess = mguess(1:4,:);
    vzguess = mguess(5,:);
    
    % Turn off display of message by fsolve
    options = optimoptions(@fsolve,'Display','off');
    
    % Solve for labor allocation and cutoff productivity levels by sector,
    % at given levels of aggregate variables
    parfor igood = 1:ngoods
        phi = vphi(igood);
        vz = mz(:,igood);
        lguess = mlguess(:,igood);
        zguess = vzguess(1,igood);
        objfun_l = @(vloglabor) residual_for_labor(Par,phi,vz,agg,vloglabor,zguess);
        l_sol = fsolve(objfun_l,lguess,options); % Solve for log of labor allocation
        mlabor(:,igood) = exp(l_sol); % Store solution of labor allocation
        objfun_z = @(logztildeH) residual_for_ztilde(Par,phi,vz,agg,exp(l_sol),logztildeH);
        z_sol = fsolve(objfun_z,zguess,options); % Solve for log of cutoff productivity
        vztildeH(igood) = exp(z_sol); % Store solution of cutoff productivity
    end
    % Compute total labor demand in each country
    LH = sum(mlabor(1:2,:),'all');
    LF = sum(mlabor(3:4,:),'all');
    % Compute trade balance as Home's trade surplus
    vstrike = (zH_expost<vztildeH); % Vector of dummies, =1 when there is strike
    moutput = mlabor.*kron([zH_expost.*(1-kappa).^vstrike;mz(2,:)],[1;1])./[1;tau;tau;1];
    sector_output = [(moutput(1,:).^((rho-1)/rho)+moutput(3,:).^((rho-1)/rho)).^(rho/(rho-1));(moutput(2,:).^((rho-1)/rho)+moutput(4,:).^((rho-1)/rho)).^(rho/(rho-1))];
    sector_price = ([agg(1);agg(2)]./sector_output).^(1/sigma);
    mprice = repmat(sector_price.*sector_output.^(1/rho),2,1)./moutput.^(1/rho)./[1;tau;tau;1];
    home_trade_surplus = sum(mprice(2,:).*moutput(2,:)-mprice(3,:).*moutput(3,:));
    res = [LH-Par.LS LF-Par.LS home_trade_surplus];
end